/*  tools.h  by K.Tsuru  */
// ver. 2.17
/**************************************************************
 bisection method
 bisectn.cpp
 It solves an equation f(x) = 0 in the region of [a, b].
**************************************************************/
#ifndef TOOLS_H
#include "tools.h"
#endif  // TOOLS_H

double bisection(double a, double b, double tolerance, double (*f)(double)) {
    if (tolerance < 0) tolerance = 0;
    if (b < a) { double c = a;  a = b;  b = c;  }

    double fa = f(a), fb = f(b);
    if (fa == 0) return a;
    if (fb == 0) return b;
    if (sameSign(fa, fb)) {
        fprintf(stderr, "f(a) * f(b) has positive sign in \"bisection()\".");  exit(-1);
    }
    for ( ; ; ) {
        double c = (a + b) / 2, fc = f(c);
        if (fc == 0) return c;
        if (c - a <= tolerance || b - c <= tolerance) break;
        if (sameSign(fc, fa)) {  a = c;  fa = fc;  }
        else                  {  b = c;  fb = fc;  }
    }
    return (a + b) / 2;
}